Classical light dispersion theory in a regular lattice 
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We study the dynamics of an infinite regular lattice of classical charged oscillators. Each individual 
oscillator is described as a point particle subject to a harmonic restoring potential, to the retarded 
electromagnetic field generated by all the other particles, and to the radiation reaction expressed 
according to the Lorentz-Dirac equation. Exact normal mode solutions, describing the propagation 
of plane electromagnetic waves through the lattice, are obtained for the complete linearized system 
of infinitely many oscillators. At variance with all the available results, our method is valid for any 
values of the frequency, or of the ratio between wavelength and lattice parameter. A remarkable 
feature is that the proper inclusion of radiation reaction in the dynamics of the individual oscillators 
does not give rise to any extinction coefficient for the global normal modes of the lattice. The 
dispersion relations resulting from our solution are numerically studied for the case of a simple cubic 
lattice. New predictions are obtained in this way about the behavior of the crystal at frequencies 
near the proper oscillation frequency of the dipoles. 

PACS numbers: 03.50.De, 41.20.Jb, 42.25.Lc, 78.20.Bh 



I. INTRODUCTION 



The classical theory of dispersion is a subject with a 
long and noble history [l], |3] . Although the main features 
of the phenomenon can be described by treating matter 
as a continuum characterized by macroscopic quantities 
such as the electric and magnetic polarizations, it is clear 
that a truly fundamental theory has to be based on a mi- 
croscopic model of matter. We shall now try to summa- 
rize some crucial aspects of the problem in an historical 
perspective, before illustrating the new features of our 
present approach. 

By treating an elementary electric dipole as an oscil- 
lator subject to a linear restoring force, it is possible to 
obtain a simple expression for the molecular polarizabil- 
ity, i.e. the complex frequency-dependent linear coeffi- 
cient which relates the microscopic dipole moment to the 
amplitude of the incident electromagnetic radiation. In 
order to correctly apply this simple model to the descrip- 
tion of the behavior of a large system of mutually inter- 
acting dipoles, one has however to consider that the field 
acting on each microscopic oscillator cannot be simply 
identified with the macroscopic electromagnetic field in 
the medium. In fact, while the latter simply represents 
the average of the microscopic field over a region much 
larger than the intermolecular spacing, the former has to 
be carefully calculated by evaluating and summing, on 
the site occupied by the considered dipole, the retarded 
fields generated by all the other dipoles of the medium. 
This "exciting" field (as we shall refer to in the follow- 
ing, although the names "effective" or "local" field have 
also been employed in the literature) was theoretically 
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estimated by Lorentz already at the end of nineteenth 
century [3|-l6| by dividing the medium into two regions 
separated by a virtual sphere surrounding the considered 
dipole. He restricted his attention to the typical situ- 
ation in which the wavelength of the macroscopic elec- 
tromagnetic field is of a larger order of magnitude than 
the average intermolecular spacing, so that one can take 
for the virtual sphere a radius intermediate between the 
two. He then argued that the influence of the portion 
of the medium lying outside the sphere can be fairly ap- 
proximated as that of a continuous distribution of electric 
dipole moment, whereas the sum of the forces exerted by 
all the dipoles situated inside the sphere can be assumed 
to vanish in most cases on the basis of symmetry consid- 
erations. A rather similar analysis, leading to equivalent 
conclusions, was also carried out by Planck [T!]. With 
these arguments one can derive the well-known Lorentz- 
Lorenz formula ,8] , relating the macroscopic dielectric 
constant of the medium to the molecular polarizability, 
and it is thus possible to deduce an approximate expres- 
sion for the dispersion relation of an array of oscillators 
in the long-wavelength regime (which generally includes 
the optical frequencies). 

A detailed microscopic theory of dispersion in a crys- 
talline solid, although with neglect of radiation reaction, 
was formulated by Ewald [13] . He considered a rectan- 
gular parallelepiped as the unit cell of the Bravais lattice, 
and his results were subsequentlygeneralized by Born to 
more general crystal structures [l|, lllj . The mathemat- 
ical methods used by these authors (one has to keep in 
mind that the theory of distributions was not yet existing 
at that time) led however to rather clumsy expressions 
for the exciting field, which could be numerically evalu- 
ated only in the limit of an infinitely large ratio between 
wavelength and lattice constant, i.e. still essentially in 
the continuum approximation. In this way the previous 
results by Lorentz and Planck were recovered for struc- 
tures with tetrahedral symmetry. Furthermore, in the 
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case of parallelepipeds of unequal edges, Ewald was able 
to perform in the same limit a numerical calculation re- 
lating the ratio between the edges to the phenomenon of 
double refraction. Finally, Ewald extended his analysis of 
the model also to the study of X-ray diffraction [12j] , but 
he made use to this purpose of other important simplifi- 
cations which are possible only in the opposite limit of a 
radiation frequency much higher that the characteristic 
frequencies of the crystal. 

Many investigations were later devoted to the applica- 
tion of quantum mechanics to the theory of light disper- 
sion, and the results of Ewald and Born were apparently 
considered to be the final word about the problem of 
the mutual interaction of a large array of classical res- 
onators. We are going to prove that, on the contrary, a 
deeper analysis reveals important properties of this fun- 
damental dynamical system which have been for many 
decades completely overlooked. 

In the present paper we shall study a system of in- 
finitely many charged particles, subject to linear restor- 
ing forces towards their equilibrium positions at the sites 
of a regular lattice, and interacting with each other 
through the retarded electromagnetic fields. Our outset 
will therefore be similar to Ewalds's, but with inclusion 
in the equations of motion of the usual "triple-dot" radia- 
tion reaction term, which corresponds to the nonrelativis- 
tic form of the Lorentz-Dirac equation |^ [l^ . The only 
approximation that we shall use is that of small oscilla- 
tions: this will allow us still to deal with a system of linear 
equations. We shall provide a general and rigorous pro- 
cedure for the calculation of the exciting field, avoiding 
to introduce at any stage of the procedure the continuum 
approximation. This will be accomplished by a method 
involving the careful subtraction of two divergent quanti- 
ties (representing respectively the total field and the field 
generated by the dipole under consideration), which ap- 
pears to be more powerful than that used by Ewald and 
Born, and presents some formal analogy with the renor- 
malization techniques of quantum field theory. Using our 
procedure we shall show that for an infinite regular lat- 
tice there exists a continuous set of normal modes which 
describe the propagation of plane electromagnetic waves. 

A remarkable result will be that the inclusion of the 
radiation damping term in the equations of motion of 
the oscillating particles, instead of giving rise to an ex- 
tinction coefficient for the wave, as is commonly believed 
according to the standard approximated treatments of 
the model, is on the contrary essential for justifying the 
presence of undamped collective waves. Such a result 
in fact constitutes an extension to the three dimensional 
case of an analogous one already obtained by two of the 
present authors for the case of a rectilinear chain of one- 
dimensional oscillators fl^. It relies upon a remarkable 
identity which was originally formulated in a different 
context by Wheeler and Feynman [l^. These authors 
deduced it from the hypothesis of the "complete ab- 
sorber" , which they introduced in order make their time- 
reversible action-at-a-distance electrodynamics compat- 



ible with the observed phenomenon of radiation reac- 
tion. For the physical system here considered we are 
going to prove in a simple and direct way that, although 
no absorption mechanism is present in the model, this 
"Wheeler-Feynman identity" actually holds as a purely 
mathematical property of the entire class of solutions on 
which we are interested. 

By numerically studying the dispersion relations for 
the crystal, as resulting from the exact solutions of the 
model, it will also be shown that completely new fea- 
tures appear for frequencies in a region about the proper 
frequency of the oscillators. In such a region the wave- 
length can in fact become as short as the lattice spacing, 
so that the approximations adopted in the previous lit- 
erature become unavailable, and only an exact solution 
can give predictions about the behavior of the system. 
It is found that inside the interval of frequencies where 
undamped wave propagation was believed to be impos- 
sible, plane waves can actually propagate, with very low 
group velocities, along certain lattice directions and for 
appropriate wave polarization. 

II. THE MODEL 

Let us consider an infinite three-dimensional simple 
Bravais lattice, that is an array of points 

Tn = niai + n2a2 + nsaa , (1) 

where n denotes the triple of relative integers (ni , n2, ns), 
and ai , a2 , aa are a set of primitive translation vectors 
for the lattice We choose the orientation of the 

ai in such a way that = ai • a2 x 33 > 0. V clearly 
represents the volume of the primitive cell. Each point rn 
is the equilibrium position of a point particle (electron) 
of mass TO and electric charge e, which is subject to an 
elastic force of the form 

Foi = -K{zn - r„) , 

where the vector Zn represents the instantaneous coordi- 
nates of the point particle. Although it is not strictly nec- 
essary for the mathematical self-consistency of the model, 
in order to reproduce in a more realistic way the situa- 
tion of solid state physics we can assume that Tn is also 
the seat of a static positive ion of charge — e. As is tra- 
ditionally the case for the classical models of dispersion, 
we shall neglect the Coulomb interaction between the ion 
and its associated electron, since classical mechanics fails 
at such short distance scales, and we shall instead iden- 
tify the proper oscillator frequency ujq = \J Kjra with a 
characteristic excitation frequency of the optical electron 
in the atomic ground state. According to quantum me- 
chanics, it could even be possible to obtain a more real- 
istic model by associating to each atom a set of fictitious 
oscillators [17| . one for each allowed quantum transition 
from the ground state of energy Eq to an excited state of 
energy En, with proper frequencies won = (£'n — Eoj/h. 
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In order to assign to the contribution of each oscihator an 
appropriate weight, one has then to make the substitu- 
tion e^/m /o„e^/m, where the "osciUator strengths" 
fon are coefficients subject to the sum rule J2n fon — 1 
(we are considering here atoms with a single optically 
active electron). An explicit calculation from first order 
perturbation theory provides 

fon = {2m/h)ujon \ {Eo \e ■ x| S„)|^ , 

X being the position operator and e the unit vector repre- 
senting the direction of oscillation. For the sake of sim- 
plicity we shall consider in the following calculations a 
single oscillator per atom, in accordance with the original 
literature on classical dispersion, although the extension 
to the case of multiple oscillators presents no conceptual 
difficulty. 

The charge and current densities associated with the 
electron-ion pair are given respectively by 

p„(x,i) = e[53(x-Zn(i))-^3(x-r„)] (2) 
j„(x,t) = eZn{t)S^X-Zn{t)), (3) 

where d denotes Dirac's delta function. The retarded po- 
tentials generated by the charge-current density (with 
= cpn) are defined as 



and putting x„ = z„ — r„, we can write the (nonrela- 
tivistic) equation of motion of the electrons as 

mxn -Kxn + een(r„ + Xn, i) 



-exn X bn(r„ + x„, t) + g^^n , (9) 



where 



Exp ret 

E-Dl-Ct 



(10) 

(11) 



represent the exciting fields. The notation used in the 
two last equations means that the summation index m 
runs over all the values in Z'^ except 0. The last term 
in Eq. © describes the radiation reaction force, accord- 
ing to the Lorentz-Dirac prescription. See Ref. [l^ for 
a comparison between the expression of the atomic po- 
larizability resulting from this classical equation in the 
dipole approximation, and the corresponding result ob- 
tained for a two-level atom in electric-dipole interaction 
with the quantized electromagnetic field. 



a;^'-*(x) 



where 

-Drot(a;) 



1 



(27r)4c 
-|x|) _ 



1 



m), (4) 



47r|x| (27r)4 J F 

is the retarded Green function 1181 and 



(5) 



III. NORMAL-MODE SOLUTIONS 

We suppose the amplitude of the oscillations to be 
small enough, so that at every stage we can neglect all 
terms of order higher than one in the Xn and their time 
derivatives of any order. It follows that Eq. ([9]) simplifies 
to 



mx„ = -ifx„ -f een(rn, t) + 



(12) 



d'^ye 



-ik-y -fi 



(6) 



is the Fourier transform of j^. We are here using the 
four-dimensional notation so that, for instance, x denotes 
the four- vector (x'^,x), and x^ = ct, k ■ x = k^x^ = 
k • x — fc°a;°, k^ — k ■ k — li.^ — (fc°)^. Summation over 
repeated indices is always implicitly understood. The 
retarded potentials satisfy the equation 

and the Lorentz gauge condition 



0. 



Introducing then the retarded fields according to the 
usual relations 



grct 



V X A"=' 



_d_ 



(7) 
(8) 



where the retarded fields, included into e„ according to 
Eq. ([TU)) . are to be calculated in the dipole approxima- 
tion, whereby each E^* depends linearly on Xm and its 
time derivatives. Of course, Eq. ((T^ actually represents 
an infinite system of coupled linear equations, since we 
have one such equation for each n ^ 1? . We shall look 
for a global solution of the form (in customary complex 
notation) 



C exp [i{K ■ Tn — Lot)] , 



(13) 



representing a plane wave with amplitude C, frequency 
Lo and wavevector k. The parameters uj and k can al- 
ways be chosen so that w > and n belongs to the first 
Brillouin zone of the crystal. We shall proceed as follows: 
for a generic motion of the form ([T^ we shall calculate 
the resulting expression for the exciting field ej,; then by 
substituting this expression into Eq. ((T2|) we shall find 
out that the equations of motion of all the particles can 
be simultaneously satisfied, provided that k and lu satisfy 
a well defined dispersion relation. 
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In the dipole approximation, that is to first order in 
C, the charge and current densities become 

Pa{x,t) = — eC • V(5^(x — Tn) exp • Tn — wi)] 
jn(x, = — jewC(5'^(x — Tn) exp • Tn — wi)] 

whence, according to the definition ^ 

Pn{k) = -27riek-C<5(/fc°-w/c)exp[i(K-k)-r„p4) 
jn(fc) = — 27riea;C(5(/c'^ — w/c) exp — k)-rn] .(15) 

Substituting these expressions into Eq. (HJ and then ap- 
plying Eq. d?]) we get 



E'r(x,t) 



eexp ■ Fm — ^t)] 

d^k (w/c)2C- (k-C)k 



(27r)3 k2-(w/c)2 
X exp [ik ■ (x - Tm)] • 



It follows 

en(i"n,i) = eL''°*(K,a;) • Cexp [i{K ■ Tn — ojt)] , 
where the second rank tensor L'^'^' has components 



(16) 



(17) 



E 



(27r)3 k2 - (w/c)2 



X exp [i{k — /t) • r„ 



and we have used the notation L''°* • C to denote the 
vector with components U^'^Cj. By substituting the ex- 
pressions (fT3|) and (fT7|) into Eq. (fT2|) , the original system 
of infinitely many coupled equations is transformed into 
the single vectorial equation 



(18) 



which admits solutions with nonvanishing C when 



det 



(k,w) 



(19) 

The above equation determines an implicit relation be- 
tween K and w, which constitutes the sought for disper- 
sion relation of the crystal. 

Once a solution of the form (IT^ has been obtained, 
it is easy to write down the expression for the macro- 
scopic quantities which can be associated to the normal 
mode. The macroscopic polarization density is in fact 
given by the continuous function of space which interpo- 
lates the microscopic displacement vectors of the individ- 
ual dipoles: 

P(x,t) = F"^eCexp[i(K-x-wt)] . (20) 

The macroscopic fields Emac, Bmac and Dmac = Emac+P 
are then given by the solutions of the usual macroscopic 



Maxwell equations inside the (nonmagnetic) medium: 



= 



cot 



cot 



(21) 
(22) 
(23) 
(24) 



From Eq. (gT]) it follows V • E,nac = -V • P. Then by 
taking the curl of Eq. (j22p and eliminating Bmac with the 
aid of Eq. ([M)) one easily obtains 



En 



P- V(V-P). 



From this, Eq. and Eq. (|22]), we conclude 

e (a;/c)2C - k{k ■ C) 



Eniac t^ 



V k2 _ (^/c)2 

X exp [i(K ■ X — ujt)] 
e {u;/c)k X C 

V k"^- (a;/c)2 

X exp [i(K ■ x — ujt)] 



(25) 
(26) 



IV. THE WHEELER-FEYNMAN IDENTITY 



Since (|19p is a complex equation, it is a priori to be ex- 
pected that, in order that a real solution for uj may exist, 
the vector k must necessarily be assigned an imaginary 
component, which represents an extinction coefficient for 
the wave. A fundamental observation can however be 
made at this point, showing that this is not actually the 
case and that Eq. ([T9|) determines cj as a real function 
of the real independent variable k. To this purpose, let 
us introduce the advanced potentials A^'^*^^ defined by 
a formula analogous to (|4]) with, in place of Diet (x) , the 
advanced Green function 



.5(x" + |x|) 
An Ixl 



1 



(2^)^ 



d^k 



fc2 



■ le. 



A completely general result about the Lorentz-Dirac 
equation asserts that the self-force, given by the expres- 
sion involving the triple time-derivative of the particle 
position, is equal to the electromagnetic force exerted 
on the particle by one half the difference between the re- 
tarded and advanced fields generated by the particle itself 
(13I [20j . Using this result in the dipole approximation, 
we have that the last term of Eq. ([T2|) can be expressed 



(27) 



67rc3 



x„ = eE(r)(r„), 



where = {Af^'"'^ ± Ai^'"^^)/2 and E± = -VA^^^ - 

{d/cdt)A^. Note that the field En~^ is a solution of 
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the homogeneous (i.e. source-free) field equation and 
therefore, at variance with En^\ it is regular at the 
particle position r^- Using Eq. (f27|) and the identity 
E^f' = eL^' + eL"\ we can rewrite Eq. ^ as 



mx„ = -Xx„ + ee(,+) (r„, t) + eE^-^ (r„, t) , (28) 
where 

e(±) - V E^^' 

/ y ^n+m 

m=iO 

m 

the summation being extended to all the values m G Z'^ 
in the last equation. We have 

A^(-)(:,) = ^Aii(-)(x) 

m 

d4fce*^-"<5(fc2)e(fc0)j^(fc) (29) 



2(27r)3c 



with 



D^-\x) EE i[Act(a;)-i?adv(a;)] 



47r 2(27r)3 



d*fce'^-"(5(r)£(F') 



and 



For a normal mode solution the above threefold series 
can be evaluated by using Eqs. (|14m5p and the relation 



^exp(^k.rO = ^^<53(k-G„ 



(30) 



where the Gm, for m g Z'^ , are the points of the recip- 
rocal lattice , defined as 

Gm = y£i3kmia.j x afc . 

Here Sijk indicates the completely antisymmetric tensor 
with £123 = 1. We obtain 

p{k) = -iek • C^^<5(fc'' - w/c) 



^,5^(k-K-G„ 



j(fc) = -iea;C^^,5(fc" -w/c) 



^,5^(k-K-G„ 



(31) 



(32) 



We see that the integrand on the r.h.s. of Eq. ([29)1 
is a singular function with support on the cone k'^ — 
|kp — (fc°)^ = 0. On the other hand, for a given real 
K, the functions p{k) and j(k) are different from zero 
on this cone only when there exists m G Z'^ such that 
Lj = c\k + Gjiil = Wm('«)- This proves that for all values 
of except those belonging to the discrete set of singular 
values {oJin{i^)}jn(zz3, the "Wheeler-Feynman identity" 







holds at any point x of spacetime. On the other hand, 
the results that will be obtained in the next section show 
immediately that, for a given k, the exciting field di- 
verges for uj — Wm('*), so that none of these frequency 
values can possibly correspond to a normal mode solu- 
tion. The Wheeler-Feynman identity is therefore estab- 
lished in complete generality for all physical solutions 
expressible as linear combinations of normal modes. As 
an immediate consequence of this identity one has that 
the last term on the r.h.s. of Eq. vanishes. Recalling 
Eq. 1221), this result can also be put in the form 



(rn,i) 



A relation physically equivalent to this one was obtained 
by Oseen (although without a rigorous mathematical 
proof) already in 1916 [2]| . 
We can then write 



e^^-'(r„,t) — eL*-^-* (k, cli) • Cexp [i(K ■ r„ — ujt)] , 



with 



L^7\-,u:) = -f (-) S. 

on \ c 



P 



(33) 

dPk {uj/c)^6ij — kikj 



(27r)3 k2 
X exp [i(k — k) ■ Tjn] 



(34) 



The symbol P J indicates the principal value of the inte- 
gral. Note that for real k and cj the function L^^-* {k, uj) is 
real, since taking the complex conjugate amounts to mak- 
ing the substitution m — )> — m in the summation index. 
From these considerations it follows that the equation of 
motion (jl8|) can be rewritten as 



-L(+)(/^,w) 
m 



(35) 



and the original complex equation ()19p is converted into 
the real equation 



det 



0. 



which determines a dispersion relation between the real 
variables k and uj. Note that this remarkable result, 
which allows for the propagation of undamped plane 
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waves in the crystal, holds just as a consequence of the 
inclusion of the Lorentz-Dirac radiation reaction term in 
the equations of motion. This does not appear surpris- 
ing, when one recalls that the expression of this term 
was determined just in order to insure global energy con- 
servation for the complete system of particles and field 

MM- 



V. THE CALCULATION OF THE EXCITING 
FIELD 

A. Outline of the procedure 

The expression for the retarded field produced by an 
oscillating dipole is well-known, and it can in fact be 
obtained by explicit calculation of the integral on the 
r.h.s. of Eq. It might seem therefore that the most 

direct way of calculating the exciting field e„ would be to 
substitute such an expression into Eq. ([TU]), as was done 
in Ref. for the one-dimensional case. It turns out 
however that in the three-dimensional case the series on 
the r.h.s. of Eq. (fTO|) docs not converge in a proper sense. 
It is possible indeed to assign to the sum an unambiguous 
meaning via the prescription 



i(rn,i) 



lim 

?;-s-0+ 



,{rn,t) exp (-rirl 



but the above formula is not convenient for numerical 
computations purposes (the sum converges slowly for 
nearly vanishing 77), and furthermore gives us no insight 
into the physical content of the results. We shall there- 
fore follow a different path, which consists in converting 
the sum over the points rn into a sum over the points 
Gm of the reciprocal lattice. After some manipulations 
we shall obtain an absolutely convergent series, which in 
typical cases can be numerically computed with little ef- 
fort. Furthermore, our final expression will provide an 
explicit expansion in powers of a/X (a = V^^^ being the 
mean lattice parameter and A = 2'k/k the wavelength), 
in such a way that the well-known result provided by the 
old theories in the long- wavelength limit will appear to be 
just the zero-order approximation of the general result. 
We start from the relation 



e„(r„,i) 



lim 

X— f Fj, 



[E-*(x,t)-E-'(x,i)] , 



(36) 



E-*(x,t)-^Er(x,t) 



where 



is the total retarded field. Using Eqs. ([T6| and (pOj) we 
obtain 

F-Vxf^ = e V (^/c)^C - [(G^ + K,)-C] (G^ + n) 

^ ' ' (G„ + k)2 - (w/c)2 - i£ 

X exp [i(Gm + n) ■ X — iojt] . (37) 



Note that, according to Eq. ([25)) . the term for m = of 
the above series is equal to the macroscopic field Emac- 
This of course corresponds to the fact that, for A S> a, 
the macroscopic field just represents the average over a 
unit lattice cell of the total microscopic retarded field. 
We observe also that r„ • Gm — 27rn • m, so that 

exp [i(Gj„ + k) ■ x] 

= exp[iK ■ Vn + i(Gm + k) ■ {X- Tn)] ■ 

Therefore, substituting Eqs. and ([H]) into Eq. 
and then comparing the resulting expression with Eq. 
(HZl), we get 

VLTMk,oj) = lim / rf3 ^p.x /''^»j-P»Pj 



P" - / - 
(5^(p-H^-q)-l 



where we have introduced the dimensionless variables 
/ = a;a/27rc, q = Ka/27r, and Hm = aGm/27r. Let 
us now shortly rewrite the r.h.s. of the above equation 
as limx^oF(x), with F{x) = / d^p e'P-^f'(p). The reg- 
ularity of the function F{x) at x = descends obviously 
from the fact that en(x, t) is regular at x = rn, as can be 
seen from the very definition ([10]) of en. Recalling that 

exp(— x2/477) 
8(7r77)3/2 



exp(-?7k2 - zk • x) 



6\x) 



we can write 



hm F(x) 



lim 

i)-i-0+ 

lim 

j)->-0+ 



d^xF{x 



d^k 



fovT]^0+, (38) 

exp(— x2/477 — iq • x) 

' 8(7777)3/2 

rf3pgi(p-q)x^(p) 
-rjlc^ — ik-x 



(27r)3 



lim 

r)-i-0+ 



d^kFik + qje-"^' 



Reintroducing the explicit expression for F, we thus ob- 
tain 

Mlfi<lJ) - - hm / d'kc,,{k,qj) 
^<53(k-Hm)-l 



(39) 



with 



Cii(k,q,/) 



{h + q,){kj + qj) - /2(5,. 
(k + q)2 - p-ie 



By isolating the term for m = of the series on the 
r.h.s. of Eq. (l39l) we get the expression 



Cii (0,q,/) 



P5, 



q2 _ j2 



Af--(q,/), (40) 
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which is related to the macroscopic field inside the crystal 
according to Eq. (1^ : 



E^ac(x, t) = eV-'M'^^'iq, f) ■ C exp ■ x - ujt)] 



We can then write 



(41) 



- hm [5,,(q,/,r;)-/,,(q,/,r;)] ,(42) 

7)— ^0 + 



is reasonable to expect for frequencies not too close to 
resonance), we have / ~ g = a/A, and it is immediate 

to check that cf^ = 0((a/A)'), Cij = O ((a/A)"*) for 
A — > oo. 



B. The lowest-order term 

Introducing yn = r^/a, we can rewrite Eq. (j30p in 
terms of dimensionlcss quantities as 



with 



m=iO 

Iio{'n,ciJ) = j d^kQj(k,q,/)exp(-ryk2) . 

The series Sij and the integral are both divergent in 
the limit 77 — >■ 0. We are going to study them separately 
and to split each of them into a divergent and a conver- 
gent part. The two divergent parts must of course cancel 
each other, as we shall check directly. We shall then be 
able to express MJJ* as the sum of a finite term and of 
an absolutely convergent series. 

By expanding the function Cij in powers of l/Zc (from 
now on k stands for |k|) we can write 



^exp(2^*k.yJ=^,S3(k_H, 



(43) 



Let us then denote with F the Fourier transform opera- 
tor, which transforms the generic function /(k) into the 
function 



[F/](x) = j d^ke^'^-/(k). 

If the function / is continuous at all points Hm, we have 
from Eq. 



JO) 



„(0 



^/(H„) = /'d3kj(k)^^3(k_H^) 

m m 

= / d3k/(k)^cxp(2^zk.yj 

n 

= 5][F/](27ry„). (44) 



where, for / = 0, . . . , 3, c-^ is a homogeneous function Applying this formula and Eq. (|38|) we obtain 
of k of degree — /, and = 0(fc~^) for k 00. The 3/2 
term of degree —3 has been multiplied hy 9{k — e) [where gxp (— nH^ 

6'(x) = 1 for a; > 0, 6'(x) = for X < 0)] in order that all ^ \V 

terms appearing in the above equation be integrable in a 
neighborhood of k = 0. We have 



.(2) 



.(3) 



1 



' \ 3/2 



(45) 



qikj + qjki 



k 



f% 



kj k T 

- 2{qikj + qjk, 
4(q-k)2 



for rj 0+, where 0{ri°°) indicates a term such that 
limj^^o+ 0{ri°°)/r]^ = for any arbitrarily large N . Let 
us then decompose cl""* as 

JO)/, X _ ^ , 



where 



1 

fc4 



i^f (k) = fc,A:, -%fcV3 



2q-k((jigj-/ %) 



+('Zjfcj + qjh 

kikj 
-4— ^q • k 



4(q.k)2 



fc2 

2(q.k)2 



is a quadratic symmetric tensor satisfying P^^ (k) 
0. Since 



fc2 



(46) 



This decomposition, besides separating the terms which 
give rise to divergent contributions to Sij and lij, leads 
naturally to an asymptotic expansion of these quantities 
for long wavelengths. In fact, if we suppose that the 
refraction index n = ck/lj = q/ f is oi order unity (as it 



is discontinuous for k = 0, we cannot directly apply Eq. 
(j44|) to evaluate its contribution to the series Sij. How- 
ever, since 

|^3ke.k.x-,.Vf (k) = -|^p(2)(x)exp(-xV477) 
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we can use Eq. p4|) to obtain 

' m5^0 m 

= ^5,,(2^y„,,7) = 0(,7-). 



It follows that 



(0) 



0(rf 



(47) 



for 77 — >■ 0^, where /3|^^^ is a finite symmetric dimension- 
less tensor, satisfying the condition ^ - fi^^^ = 0. The 

numerical values of (3^^ in general depend on the par- 
ticular crystal structure considered. In Appendix A we 
derive the following equivalent expression in terms of a 
sum over the points of the direct lattice: 



.(2) 

4^--o+— |y: 



3 \ - P-t] (yn) , 2n 

i:; 15 exp(-ry„ 



(48) 



From Eqs. and (|T7)) we conclude 

m=iO 



"'3 

3 



3/2 



+ ^+0(^7°°)- (49) 

On the other hand, it is almost immediate to see that 

J d^kbf{k,r]) =0. (50) 

In fact, for obvious symmetry reasons, the l.h.s. must be 
of the form f{ri)5ij, but then the condition X^j^ii''' — 
implies /(?]) — 0. It follows that 



if^ ^ J d^kcfik)eM-vk') 



3/2 



^ J d'k eM-vk') = f ) (51) 



and so 



lim ( S, 



(0) 



r)-!-0+ 

We thus conclude that at order zero in a/ A the tensor 
M-j*- is given by 

(2/2 + q2)5„./3-g,gj 



(q,/) = (q,/) + ^-/3j) 



P'lf- (52) 



For an isotropic crystal the optical behavior at long 
wavelengths must be invariant under spatial rotations. 



This means that the symmetric tensor /3|^^ must be a 
multiple of the identity matrix, but since its trace ^ ^ /Jj-^"'' 



,(0) 



0. In such 



vanishes, we see that isotropy implies (31^ 
a case, recalling Eqs. ([20)1 and (HTj) . we derive from Eq. 
([52]) that at this order of approximation 

e„(r„ + -P(r„,i). (53) 

The above expression for the exciting field is the same 
that was obtained by Lorentz and Planck with the argu- 
ment of the virtual sphere mentioned in the first section, 
and was later confirmed through more rigorous analysis 
by Ewald and Born. 

The same matrix M^^^ of Eq. ((5^ of course also 
provides the zeroth order approximation for the tensor 
M(+) = Vi'^+\ Therefore recaUing Eq. ^ one obtains 



En 



-P-/3(") -P, 



(54) 



where uip = e/ \/raV is the so-called "plasma frequency" 
of the material. The above equation can be put into the 
form P = x(w) ■ Einac, where the tensor 



x(^) 



-w2_^2/3 



(55) 



represents the electric susceptibility. The dielectric func- 
tion can then by obtained as £(w) = l-l-x(w). It appears 
from Eq. (I55|) that at the present order of approximation 
the lattice of resonators behaves in general as a biaxial 
crystal. Let us denote with /3|°'' (i — 1,2,3) the three 
(in general distinct) eigenvalues of the symmetric tensor 
^(0), with 



£,(a;) = 1 + 



(0) 



the eigenvalues of £(w), and with ni{uj) — ^ £i{io) the 
indexes of refraction for transversal waves with electric 
fields polarized along the mutually orthogonal directions 
of the corresponding eigenvectors. We then have that, 
for i ^ j, the (approximately) frequency independent 
quantities 



1 



1 



1 



1 



that were introduced by Havelock [22| as a measure of 
the phenomenon of structural double refraction, can be 
calculated according to our model as 



D^J - A 



(0) 



(56) 



According to Eqs. (|T7)) or P5|) these parameters depend 
solely on the lattice structure. 
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As we mentioned in Section I, Ewald was able to de- 
vise a method to calculate when the unit cell of the 
lattice is a rectangular parallelepiped Although his 
formulas look considerably more complicated, we have 
checked that the numerical results, that he obtained for 
particular values of the ratio between the edges of the 
cell, are in very good agreement (considering the tools 
available at that time for numerical computation) with 
the general result provided by Eq. (|56)) . 

In the following subsection we are going to show that 
Eq. ([5^ has to be significantly corrected when the size 
of the lattice parameter is not negligible with respect to 
the wavelength. It will be found that it is no longer 
possible in that case to give for the dielectric tensor a 
general expression as a function of the frequency only. 



since, as was already recognized by Ewald [T^, the form 
of the dispersion relation involves both the direction of 
polarization and the direction of the wavevector k. Our 
analysis will allow us to study quantitatively in detail the 
behavior of such a dispersion relation. 



C. The complete solution 

Since c-j\ as well as c-^\ are odd functions of k, they 
clearly give no contributions to either Sij or Jy . Let 

(2) 

US now consider c,-,- . If we introduce the fourth order 
completely symmetric tensor 



fc2 
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Hjhii^) = hkjkhh ~ '^{Sijkhh + Sihkjki + Sukjkh + Sjhhki + Sjihkh + Shihkj) 



satisfying the condition J2i Puiii (k) = 0, we obtain after some simple algebraic manipulation 



c!-Hk,q,/) 



1 



15fc2 
4 



[3(7.9, -%(q' + 10/2)] 
\5,,Pl^\^)quqi - Qqn (^.^'^(k) + q,P^l\\4) + (k)(7/2 - Sq^)' 



7fc4 

'J^Ptjliii'^)9hqi- 



By integrating Eq. (|45p with respect to 77 we obtain 



E^^£^~- + -+o(,n. (57) 

m#0 ™ 



whence, using Eq. 



where a is a finite dimensionless constant dependent on 
the crystal structure. Similarly, if we put 



drf 



^6(|,(H^,,7) - 5:4^\(H^)exp(-^HL) 



^?^(k) 



It follows that 



integration of Eq. (H71) provides 



J2 ^ (H™, 77) ^ ^ - ^^if + 0(j^n , (58) 



where the integration constants 7|j2; are completely sym- 

where (3)^' is another finite symmetric tensor such that , • , ^- 1^ ■ ^^ (n) n ^ 

metric tensors satisfymg the condition i^i^l^i^i = tor 

Li Pa = 0- Let us then define n = 0, ±2. We show in Appendix A that 

(4) 

'Pexp(-ry^) (60) 



(59) 



,(2) 



We have 



/ 



^3/2 



(0) _ 2^ 1- \ ^ Pjjhl iyn) ^1 _ 2n 



TT72 4^^z(^)exp(-xV4^) 



I677I1/2 y 



(_2) _ 105 



P^^^ fv ) 



(61) 



n#0 



10 



A preliminary determination of "f^jj^i and J^jj^j^ with the 
aid of the two above formulas can considerably improve 

(2) 

the accuracy in the numerical calculation of jljl^i accord- 



ing to Eq. (HH). From Eqs. §7}, ([58]) and §9!i it follows 
that 

Sl? ^ ^ c|f(H^,q,/)exp(-77HL) 

2 7r3/2 

[3q,q, -S.,ie + 10. f)] 



with 

Mif(q,/) 



15 ^ 
-Ml^\qJ) + Oiv), 



H?M<1'^<11 ■ (62) 
By the same argument used to deduce Eq. ([50]) we have 

J d^kbfik,7^)=0. (63) 

In a similar way one finds that 

J d^kb'^lik,r,)=0. (64) 

In fact, since the l.h.s. is an invariant tensor, it must be 
of the form f{r]){6ij6hi + SihSji + SaSjh), but then again 
the condition ^iiM ^ implies f{ri) = 0. It follows 
that 



I d3kcg'(k,q,/)cxp(-,7fc2) 



= [3q^qj-S,j{q^ + 10f)] / d^k 



2^3/2 



exp(— T^fc^) 



15/c2 



[3g,gj-%(q' + 10/')] , 



so that 



lim I S, 

r)-!-0+ 

Finally, we have 



lim V Cjj (Hm, q, /) exp {-rjUl 

= J2 c.,(H^,q,/) = 5.,(q,/), 

lim / d^kcij(k,q, /)exp(-77fc2) 
= / cfk c^j (k, q, /) = I^j (q, /) , 



since both the sum and the integral on the r.h.s. of the 
two above equations are absolutely convergent. An ex- 
plicit calculation carried out in Appendix B shows that 



Jy- (q, /) - -t-TT^fS^j = a4 ) (q, /) , (65) 



4 

in accordance with Eq. (|33]) . We can therefore conclude 
that 



M-*(q, /) = M,(+)(q, /) - z-7rV'% , (66) 



with 

m(+) (q, /) = m(°) (q, /) + Af (q, /) - S., (q, /) . (67) 

As we had anticipated in the previous section, the ma- 
trix m[^^ is symmetric and real for real q and /, whereas 

Af|- ' just cancels the radiation reaction term in the equa- 



tion of motion. Since m[^^ = VL)J' , we can rewrite Eq. 



-(+) 



(155]) as 



MW(/«a/27r,a;a/27rc) • C = 



(jJn — UJ 



or, in dimensionless variables, 
r(+) 



— M^+^(q,/).C-(/o^-/^)C, 



(68) 



where /o = a;oa/27rc, and rc = e'^ / ATrmc^ is the classical 
electron radius. 



VI. THE SIMPLE CUBIC LATTICE 

In order to apply the general theory developed un- 
til now to a concrete situation, let us consider in more 
detail the particular case of a simple cubic lattice, for 
which the primitive translation vectors of Eq. ([1]) are 
given by — aUi for i = 1,2,3, being the unit 
vectors of the three coordinate axes and a the lattice 
parameter. For this lattice one has Hm — m, and 
the first Brillouin zone corresponds to the cubic region 
— 1/2 < qi < 1/2. Obvious symmetry considerations im- 
ply that = iS^f = p^'^^Sij. Then from the 

general condition J2i PiP — Si Pu^ = one deduces 
that /3(o) = = 0, and so 

The symmetry also implies that 7^"^^ = X^^^'Eijhi — 

I'^'^^'Dijhu where D^jhi = kjShi + ^ihSji + &iiSjh and 
the non-tensor Eijhi is defined so that Eijhi — 1 for 
i ~ j ^ h = I, Eijhi ~ otherwise. Since Dnhi — 



hi 



/(») 



5m, the condition 'Zilzlhi 



implies 



57*^"'. Equation (|62]) thus becomes 



M^fiqJ) = ^[3q,q,~6,M' + W)] 

+47(2)(29,9,+q25y--5g2*y). 
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A numerical calculation based on Eq. ([57)) shows that while from Eqs. (|59II61I) one obtains 

2^3/2 



a = lim 



E 



n\= 8.913633, 



J 



v(-2) - L,(-2) _ _21_ 



2^"" " 32^3 



V - 0.0751838 



v(o) - 



v(2) 



^ lim y^^e— ^= 0.185800 



4t f E 'i'li"'' ^"""^ 
10 ,^0. 1^^^ |m|6 



+ -y(o)^_^(-2)^2y2 9^^ 0.2780310. 



For application in numerical computations, the above results can also be put in the useful form 



<^(q,/) = - E 

|mt<L 



(TOj +g,)(mj- - /2^^j- ^ 

(m + q)2-/2 + ^^^^T 



+ (q2 - 2/2)%/3 - 2gf5,,] - 2C(i)(2g,(7, + ~ hql^,,) 

A-TT 

+ Yo5l ["9«*(3q' - 7/') + (|q|' + 2Wf + 70/^)%] + 0(1-') , 



r 



where we have introduced the functions 



A{L) 
B{L) 
C{L) 



E 1 



|m|<L 



niT^O, |m|<L 



mT^O, |m|<L ' ' 



with /3 = a/5 - 27(2) 1.226665. 

As an example of application of these formulas, we 
have numerically calculated a few dispersion curves for 
an ideal simple cubic lattice with parameter a = 2rB 
and proper frequency of the oscillators luq — lub, where 
rg = Airh^/me^ and = {e^ /AirYm/h^ are the Bohr 
radius and the Bohr frequency respectively. In Fig. 1 
we show the dependence of |q| — Ka/ii: as a function 
of w/wo, as determined according to Eq. (|55| for three 
different directions of propagation of the plane waves. 
The results can be compared with those obtained from 
the standard formula valid in the long wavelength limit, 
which, as can be seen from Eqs. and ([5^. is given 
by 



(2/2 + q2)C/3-q(q-C) 



ira 



P 



One can see that, according to the above equation, C 
has necessarily to be either orthogonal or parallel to q. 



In the first case (transversal polarization) one obtains the 
dispersion relation 



1/2 



- Cj2 - cj2/3 



while the second case (longitudinal polarization) can oc- 
cur only when lu = lui = 



lIq + (2/3)cl;2^ independently 

of the value of k. These relations are represented in the 
figure by the thick solid lines. The curve corresponding 
to transversal polarization presents a vertical asymptote 

for iO ~ UJa = 



—u}p/3, whereas that for longitudi- 
nal polarization is obviously just a vertical straight line at 
Lo = uji- According to the standard theory no propagation 
is possible in the frequency interval between uja and w;. 
A further dispersion curve for transversal polarization is 
present for uj > uii, according to both the approximated 
and the exact theory, but its onset in the figure is prac- 
tically indistinguishable from the frequency axis. Note 
that, for the particular numerical values that we have 
considered, one has rc/ira — (wpa/27rc)2 ^ l/(27r x 137^) 
and /o ^ 1/{tt x 137), so that ujp/ujo = ^7r/2 = 1.253, 
LOa/iOo = ■\/l — 7r/6 ^ 0.690 and wi/uja = ^1 + 7r/3 = 
1.431. 

As one could expect, one can notice that the predic- 
tions of our exact theory depart in a significant way from 
those of the approximated one as soon as naj^-K becomes 
of order of magnitude comparable to unity. This situa- 
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FIG. 1: Dispersion relation for three different directions of propagation of a plane wave in a cubic lattice with a = 2rs and 
cjo = u)B- The polarization for each branch of these curves is described in the text. 



tion of course corresponds to the fact that the lattice pa- 
rameter is no longer negligible with respect to the wave- 
length. A new feature which is revealed by the exact 
theory is the appearance, even for a perfect cubic lattice, 
of an anisotropic behavior, consisting in a remarkable 
dependence of the dispersion relation on the direction of 
propagation and polarization of the plane waves. The two 
thin solid curves at the opposite sides of the figure, which 
both refer to the direction of propagation [100], corre- 
spond to transversal (the lower frequency one) and lon- 
gitudinal polarization respectively. If one considers the 
dependence of w as a function of k which is described by 
the former curve, one can see that the frequency reaches 
a maximum — which corresponds to a zero of the group 
velocity duj /dn — and then slowly decreases over a wide 
region of the k axis extending until the edge of the Bril- 
louin zone. This means that there exists a relatively small 
interval of the frequency axis with the property that, for 
each Lo belonging to this interval, there are two transver- 
sal modes with direction of propagation along a crystal 
axis and different values of the wavelength. 

A similar phenomenon occurs also for the lowest- 
frequency branch of the dashed curve, which corresponds 
to a wave again polarized along a crystal axis, but prop- 
agating along the [110] direction. A completely different 
dispersion relation is shown instead by a wave propa- 
gating in the same direction, but transversely polarized 
along the [110] direction. This is illustrated by the in- 
termediate branch of the dashed line. For this curve the 



frequency monotonically increases with n over the en- 
tire Brillouin zone, but with a slope — and therefore a 
group velocity of the associated waves — which drops to 
a relatively small value when ui approaches and then ex- 
ceeds uja- At the edge of the Brillouin zone this branch 
joins continuously, at a point with vertical tangent, the 
third branch of the same curve, which corresponds to 
longitudinal polarization. This branch extends up to w;, 
where a new transversal branch begins. It follows that, 
at distinction with the modes propagating along the [100] 
direction, there is no forbidden frequency zone for waves 
propagating along the [110] direction. Of course, no hint 
about the existence of such a phenomenon could have 
been derived from the standard theory. A very similar 
behavior is shown by the dash-dot curve, correspond- 
ing to waves propagating along the [111] direction. Here 
again the branch with transversal and that with longi- 
tudinal polarization join continuously at the edge of the 
Brillouin zone. Note that, whereas a wave propagating 
along the [110] direction can be transversely polarized 
along two directions (the [001] and [110]) with distinct 
symmetry properties and so also with different disper- 
sion relations, the two possible transversal polarizations 
are degenerate for waves propagating along either the 
[100] or [111] directions. 

For all the curves displayed in Fig. 1, the polarization 
vector C is either parallel or orthogonal to the wavevec- 
tor q. We would like to point out however that this fact 
is not a general property of the the solutions of Eq. (|68l) , 
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but is rather a consequence of the special symmetry of 
the particular directions of propagation that we have con- 
sidered. Note also that the continuous links between the 
transversal and the longitudinal branches of the curves 
can be seen as an obvious consequence of the particular 
degeneracy of the matrix M(+^(q, /) when q = (i, i, 0) 
and q = (i i i). 



VII. CONCLUDING REMARKS 

In this paper we have presented a new refined version 
of the classical dispersion theory of electromagnetic ra- 
diation in a crystalline solid, which has been completely 
described at a microscopic level as an infinite regular ar- 
ray of charged oscillators, without ever making use of the 
continuum approximation. It has been shown that, when 
the wavelength is comparable with the lattice parameter, 
the predictions provided by our calculations differ in a re- 
markable way from those derived from the old standard 
approximated formulas. In particular, in the region near 
the proper frequency of the oscillators, which was be- 
lieved to represent a forbidden gap for light propagation, 
the phenomenon of "slow light" (i.e. light propagation 
with small group velocity) is found instead to take place 
along appropriate crystal directions. 

The model we have studied is however interesting in 
itself also from a purely theoretical point of view, inde- 
pendently of its phenomenological implications. It is in 
fact an exceptional case of a completely solvable model 
(although in the dipole approximation) in which the radi- 
ation reaction force acting on classical charged point par- 
ticles is fully taken into account. Actually, it turns out 
that the inclusion of this force has fundamental implica- 
tions on the qualitative properties of the solutions, allow- 
ing for the existence of undamped collective oscillations. 
Another interesting feature is that, although the model 
was formulated in terms of the usual Maxwell-Lorentz 
electrodynamics with Lorentz-Dirac selfinteraction, it is 
also compatible with the action-at-a-distance electrody- 
namics which was proposed by Wheeler and Feynman in 
15|. In fact, when the Wheeler-Feynman identity 



1945 

is verified, the two theories lead to the same equations of 
motion for the charged particles. 

Of course the investigations that we have presented 
here can be developed and extended in several directions. 
One of these, which appears to be particularly signifi- 
cant also from a phenomenological point of view, is the 
study of the behavior of a semi-infinite lattice occupying 
only one half of the full three-dimensional space. This 
investigation should clarify the relationship between the 
normal modes of the crystal — which we have described 
in the present work — and the observable electromag- 
netic radiation propagating in the free half-space, i.e. the 
phenomena of refraction and reflection at the surface of 
the crystal. In this way it is to be expected that classical 
fundamental results, such as the Ewald-Oseen extinction 
theorem (23l - [26| . will be justified on the basis of the de- 



tailed microscopic dynamics of the complete system. 



Appendix A: Series evaluation 

Let us introduce the operator Cr, which acts on the 
generic function /(k) by operating the convolution with 
a gaussian of width y/2T: 



[a/](p) = / 



,3u exp(-|p~k|V4T) 
8(7rr)3/2 



/(k). 



According to Eq. ([55)1 . if / is continuous in p, then 
\miT-^Q+[Cr f]{p) — /(p)- Furthermore we have for any 
/ 

[Fa/](x)=e-^-'[P/](x). 

The function 6|°''(k, 77), defined by Eq. (|46l) . is contin- 
uous for k ^ and we have 



8(7rr)3/2 ^3 



C(k,^)=0, 



Therefore, applying Eq. (|44l) to the function Crbl^\ we 



can write 



T->0 + 

Using the formula 
1 



= lim V[a&g^](H^,?7) 

m 

= hm V[Pa6lf](27ry„,77) 

n 

= hm Ve-^^"[F6lf](2^y„,77). 



(Al) 



d^k- 



1 



(27r)3 J k2 Attx ' 

with standard arguments it is easy to show that 

^ [F6(")](x,,?)^--i^/^f (x)+0(7y-) 



for ?7 — > 0+, x 7^ 0, while according to Eq. ([50)) 

Therefore from Eq. (jAip one deduces immediately Eqs. 
(|T7l) and ([15)) . In a very similar way, by the same argu- 
ment we used to deduce Eq. ([64)) we have 



[Fb%]iO,v) = 0, 
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whereas one can show that 

1 .^,(0),. . 15 



(2^)3 "ijhl 



for ry — !• 0+, X 7^ 0. Applying again Eq. (|44p we then 
have 



r— ^0+ ■'^ — ^ 
n 

r— ^0+ — ' 
n 

= 7lS^-7Uf'? + 0(r), (A2) 

with jI'^Ii and 7-j7i';'' given by Eqs. ([60]) and (|6T|) respec- 
tively. Hence, by integrating with respect to 77, Eq. ([5^ 
is finally obtained. 

Appendix B: Calculation of lij 

We have 

4(q./) = / d^kcy(k,q,/) 

= lim I / d3kc,j(k,q,/) - ^L^J,^ 



47r^ 
5" 



(Bl) 



Let us put q — |q|. Since Cij is a second-rank tensor, 
from invariance considerations it is readily found that we 
must have 

/ d^kc,,(k,q,/) = lAiqJ)(^~5,,) 
Jk<L ^ \ 1 / 



with 



k<L 



d^kc,j{k,q,f) ^2 



A{qJ) : 

Jk<L 



g2 



After performing a rotation on the integration space k, in 
such a way that the unit vector U3 of the third coordinate 
axis is directed along the vector q, we obtain 



A{qJ) 



- f 2 



— P — is 



-0(L-|k-gU3| 



= -fLh{q/L,f/L) 
+L'h{q/LJ/L) 



(B3) 



where 



-TrL'-2fLhiq/LJ/L), (B4) 



k"^ — P — ie 



We have 

5(0,/) = 3A(0,/) 
4 



= -ttL-" - Stt/^ \ L + dk 



1 



3 



(B5) 



for L 00, where we have used the relation 



1 1 f+°° 1 

dk—: - dk- 



k'^ — p — ie 2 k"^ — p — ie 

+0{L-') ^ ^ + 0{L-') 



We have also for q = 0{L-^), f = 0(1-^) 



d 



h{qj) = / d-^k 



(5(1 - |k - gual) fca - g 

fc2„/2 |k-gU3| 



d-'k 



k"^ + 2qk3 + q"^ - P k 



2tt / d^ 



1 + 2ge + 9' - 
= -y.-+0(X-), 

whence 

h{q/L, f/L) = h (0, f/L) + 0(L-2) . (B6) 
In a similar way we obtain 

J-i l + 2qi + q^-p 



IQtt 



g + 0(L-3)^ 



so that 



57r g 



/2(g/L, f/L) = /2(0, //L) + + 0{L-^) . (B7) 
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From Eqs. (|B3IIB7I) we get 



A-jr 

B{qJ) = — -87: fL^iA7:'f + 0{L-'), 



A{qJ) 



9 3 V 5 



then by substituting these expressions into Eqs. (|B2|) and 
(|BT|) we finally obtain Eq. ((^ . 
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